import numpy as np
from scipy.optimize import minimize
import emcee

def likelihood(params,t,y,e):
    n=len(t)
    y.resize((n,1))
    y=y#-np.mean(y)
    C=np.empty((n,n))
    L=np.full((n,1),1)

    for i in range(n):
        for j in range(n):
            if i==j:
                C[i,j]=C[i,j]+e[i]**2
            C[i,j]=C[i,j]+params[1]**2*np.exp(-abs(t[i]-t[j])/params[0])

    C=np.mat(C)
    L=np.mat(L)
    Cperp_in=C.I-C.I@L@(L.T@C.I@L).I@L.T@C.I
    #p=np.linalg.det(C)**(-1/2)*np.linalg.det(L.T@C.I@L)**(-1/2)*np.exp(-y.T@Cperp_in@y/2)
    #prior=-1/2*np.log(params[0])-np.log(params[1])
    #return 1/2*np.log(np.linalg.det(C))+1/2*np.log(np.linalg.det(L.T@C.I@L))+y.T@C.I@y/2-prior
    #return 1/2*np.log(np.linalg.det(C))+y.T@C.I@y/2
    return np.linalg.det(L.T@C.I@L)**(-1/2)

def PRH(t,y,e,method='MAP'):
    initial_params = [500,0.14]
    soln=minimize(likelihood, initial_params,method="L-BFGS-B",args=(t,y,e),
                  bounds=[(1,1000),(0.01,1)])
    return soln.x

def DRW_process(t,tau,SF,m):
    '''Return a DRW light curve

    parameters
    ----------
    t  : array_like
         epcohs
    tau: float
         input timescale
    SF : float
         input StructureFunction_infty(=sqrt(2)*amplitude)
    m  : float
         mean of the ligh curve

    Returns
    -------
    DRW_process:ndarray
                a simulated light curve
    '''   
    r=np.diff(t)/tau
    ls=[np.random.normal(m,SF/1.414,1)[0]]
    for i in range(len(t)-1):
        if r[i]<0:print('Error:时间序列未排序');return 
        stdev=(1-np.exp(-2*r[i]))**0.5*SF/1.414
        loc=ls[i]*np.exp(-r[i])+m*(1-np.exp(-r[i]))
        ls.append(np.random.normal(loc,stdev,1)[0])
    return np.array(ls)

t=np.sort(np.random.uniform(0,2922,445))
y=DRW_process(t,292,0.2,18)
e=np.full_like(y,0.03)
s=np.random.normal(y,e)
#save=np.c_[t,s,e]
#np.savetxt('lc.txt',save,fmt='%.4f')
#print(PRH(t,s,e))
'''
initial=[500,0.14]
ndim, nwalkers=len(initial),16
sampler=emcee.EnsembleSampler(nwalkers,ndim,likelihood,args=(t,y,e))
#Running burn-in
p0=initial+1e-4*np.random.randn(nwalkers,ndim)
p0,lp,_=sampler.run_mcmc(p0,125)    
#Running production
sampler.reset()
sampler.run_mcmc(p0,500)
#Markov chains
t_chain=np.exp(np.sort(-sampler.flatchain[:,0]))
s_chain=np.exp(np.sort(sampler.flatchain[:,1]))

plt.hist(t_chain)
plt.show()
'''
for i in np.linspace(10,1000,20):
    for j in np.linspace(0.01,0.7,20):
        print(likelihood([i,j],t,s,e))
